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ABSTRACT 


This  paper  describes  the  computer  implementation  of  a  two-and-one-half 
dimensional  (2.5D)  constant  density  prestack  inversion  formalism  with  laterally  and 
depth-dependent  background  propagation  speed.  This  is  a  Kirchhoff-type  inversion, 
summing  a  line  of  receiver  data  over  traveltime  curves  in  the  depth-dependent 
background  medium  with  weights  determined  from  the  Bom/Kirchhoff  inversion 
theory.  This  theory  predicts  that  the  output  will  be  a  reflector  map  with  peak 
amplitudes  on  each  reflector  being  in  known  proportion  to  the  angularly  dependent 
geometrical  optics  reflection  coefficient.  The  2.5D  feature  provides  for  out-of-plane 
spreading  correction  consistent  with  the  prescribed  background  medium.  The  method 
is  applied  to  a  synthetic  data  set  and  to  an  experimental  data  set  generated  at  the 
Seismic  Acoustic  Laboratory  at  the  University  of  Houston  under  support  of  the 
Marathon  Oil  Company.  The  graphical  output  demonstrates  the  validity  of  the 
formalism  as  a  Kirchhoff  migration.  Parameter  estimation  for  the  experimental  data 
was  less  successful,  partially  due  to  problems  with  amplitude  control  in  the  original 
experiment  and  partially  due  to  the  limited  aperture  of  the  common  shot  data,  thereby 
suggesting  that  a  common  offset  inversion  might  be  more  useful  for  parameter 


INTRODUCTION 


The  purpose  of  this  paper  is  to  describe  a  computer  implementation  of  a  two- 
and-one-half  (2.5D)  dimensional  constant  density  prestack  inversion  formalism  with 
laterally  and  depth  dependent  background  propagation  speed.  We  refer  to  this 
inversion  procedure  as  “Kirchhoff  inversion”  because  of  the  striking  similarity 
between  this  method  and  Kirchhoff  migration  in  a  variable  background  medium.  The 
major  difference  is  a  spatial  weighting  derived  from  the  inversion  theory. 

In  this  implementation,  this  medium  is  made  up  of  constant  velocity  layers  with 
curved  interfaces,  although  this  specialization  is  not  a  requirement  of  the  theory.  The 
origins  of  this  approach  to  seismic  inversion  can  be  found  in  the  work  of  Cohen  and 
Bleistein  (1979),  Bleistein  and  Gray  (1985),  and  Cohen  and  Hagin  (1985).  Beylkin 
(1985)  proposed  a  method  of  operator  inversion  that  generalizes  the  results  of  those 
papers.  An  extension  of  that  theory,  consistent  with  the  earlier  work,  was  developed 
in  the  papers  of  Cohen,  Hagin  and  Bleistein  (1986);  and  Bleistein,  Cohen  and  Hagin 
(1987).  The  extension  revises  the  Beylkin  inversion  operator  in  two  ways.  The  first 
modification  provides  a  reflector  map  instead  of  a  velocity  model,  with  the  peak 
amplitude  on  each  reflector  being  proportional  to  the  geometrical  optics  reflection 
coefficient.  This  work  has  its  origins  in  a  series  of  papers:  Bojarski  (1967,  1968), 
Mager  and  Bleistein  (1978),  and  Cohen  and  Bleistein  (1979).  This  modification 
provides  a  quantification  for  the  more  abstract  construct  of  a  “wavefront  set”  of  the 
pseudo-differential  operator  approach  to  inverse  scattering  following  Beylkin  (1985). 
This  modification  is  a  matter  of  necessity  for  seismic  data,  which  tends  to  be  high 
frequency  data  for  the  length  scales  of  the  seismic  experiment.  From  such  band- 


limited  data,  it  is  not  practical  to  extract  information  about  the  velocity  (medium 
parameter)  field  itself,  but  only  about  the  “discontinuity  surfaces’’  —  the  reflectors  — 
of  those  parameters. 

The  second  modification  provides  a  2.5D  inversion.  That  is,  it  allows  for  the 
processing  of  a  line  of  data  to  produce  a  two  dimensional  reflector  map  with 
amplitudes  that  approximate  the  effect  of  the  out-of-plane  spreading  of  the  re  jpO*»SC  to 
a  three  dimensional  point  source.  The  method  assumes  that  the  subsurface  has  two 
dimensional  variation  only,  with  the  data  line  being  a  dip  line  of  the  subsurface. 

The  inversion  operator  presented  in  these  papers  is  based  on  the  Bom 
approximation  for  the  wave  upward  scattered  from  the  inhomogeneities  in  the  earth. 
However,  an  analytical  proof  of  the  validity  of  the  inversion  formalism  as  applied  to 
Kirchhoff-approximate  data  has  been  presented  in  Bleistein  (1987a,  1987b,  1989). 
This  proof  partially  overcomes  the  “small  perturbation”  constraint  of  the  motivating 
Bom  approximation.  The  proof  shows  that  if  the  background  medium  above  a 
reflector  is  known  (or  the  background  is  close  to  the  true  medium  above  the  reflector) 
then  the  reflector  will  be  properly  positioned  (approximately,  for  a  close  background). 
Furthermore,  the  peak  output  on  the  reflector  is  linear  in  the  angularly  dependent 
geometrical  optics  reflection  coefficient,  from  which  it  follows  that  the  change  in  earth 
parameters  across  the  reflector  being  imaged  need  not  be  small. 

As  an  alternative  to  the  use  of  the  Bom  approximation,  one  can  start  with  the 
Kirchhoff  approximation  for  a  single  reflector  and  develop  an  inversion  operator  based 
on  this  model.  This  method  leads  to  the  same  inversion  formalism,  as  might  be 
expected  by  the  method  of  proof  described  above.  An  implementation  of  this 
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approach  for  common  offset  inversion  is  presented  in  Sullivan  and  Cohen  (1987). 
Docherty  (1987a,  1987b)  also  applied  this  approach  to  2.5D  zero  offset  inversion  in  a 
medium  comprised  of  many  layers.  The  first  author  of  this  paper  used  the  same 
approach  to  confirm  the  inversion  for  prestack  common  shot  2.5D  inversion  (Dong, 
1989,  Dong  and  Bleistein,  1989)  and  to  verify  equation  (50)  in  Bleistein,  Cohen  and 
Hagin  (1987)  for  this  case.  The  computer  implementation  was  then  developed  using 
Docherty’s  zero-offset  inversion  code  as  a  point  of  departure.  For  details  of  the 
derivation  of  the  inversion  operator,  the  reader  is  referred  to  those  references. 

The  computer  implementation  was  tested  on  ray-theoretic  data  synthetically 
generated  with  Docherty 's  (1987a)  CSHOT  program.  In  addition,  tests  were  carried 
out  on  physical  model  data.  Physical  model  data  are  useful  for  testing  and  comparing 
seismic  data  imaging  techniques  (migration  or  inversion).  The  physical  model  data  are 
actual  recorded  wavefields  and  contain  all  wave  effects  including  lateral  waves,  near 
field  effects,  mode  conversions,  and  diffractions.  Seismic  data  modeling  and  imaging 
methods  are  based  on  a  theory  that  incorporates  simplifying  assumptions  about  the 
wavefield.  If  an  imaging  procedure  is  based  on  the  same  theory  as  the  modeling 
procedure,  the  imaging  procedure  is  merely  the  inverse  of  the  modeling  procedure. 
While  this  is  an  advisable  first  test  on  an  inversion  formalism,  the  imaging  might  work 
perfectly  on  synthetic  data  from  the  modeling,  but  might  not  work  on  field  data. 
Beyond  numerically  generated  data,  physical  model  data  provides  the  next  level  of 
test;  the  data  are  real  wavefields.  Since  the  models  are  simpler  than  the  real  earth  and 
th^  physical  parameters  are  known,  physical  model  data  can  be  used  to  verify  imaging 
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techniques. 


Our  physical  model  data  was  generated  at  the  Seismic  Acoustic  Laboratory  at  the 
University  of  Houston  under  support  of  Marathon  Oil  Company,  who  uses  such 
physical  model  data  to  evaluate  imaging  (migration)  processing  by  contractors.  Since 
the  result  is  known  beforehand,  this  can  be  input  as  the  migration  background  velocity 
function.  The  migration  techniques  can  be  compared  independent  of  the  velocity 
analysis.  Marathon  donated  a  physical  model  data  set  to  the  Center  for  Wave 
Phenomena  so  that  we  might  try  our  inversion  on  the  data.  The  model  is  structurally 
complicated  enough  to  warrant  prestack  inversion.  Furthermore,  the  model  was 
sufficiently  two  dimensional  to  make  the  application  of  Dong’s  code  to  this  data 
practical.  That  application  was  carried  out  by  the  second  author  of  this  paper.  An 
extensive  discussion  of  the  analysis  of  this  data,  including  a  search  for  mode  conver  d 
waves  and  lateral  waves,  can  be  found  in  Emanuel  (1989a,  1989b). 

Parameter  estimation  for  the  synthetic  data  set  showed  errors  of  about  7%. 
Parameter  estimation  for  the  experimental  data  was  less  successful,  partially  due  to 
problems  with  amplitude  control  in  the  original  experiment  Both  data  sets  suffered 
.Tom  the  limited  aperture  of  the  common  shot  data,  which  degrades  amplitude 
accuracy,  thereby  suggesting  that  a  common  offset  inversion  might  be  more  useful  for 
parameter  estimation.  Tests  with  a  common  offset  inversion  program  are  now  in 
progress  and  will  be  reported  in  a  future  paper.  Earlier  tests  on  wide  aperture 
synthetic  data  with  fine  sampling  produced  satisfactory  amplitude  estimates  and 
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convinced  us  that  the  theory  is  valid. 


THE  CXZCS  ALGORITHM 

For  a  2.5D  acoustic  model  comprised  of  constant  velocity  layers,  separated  by 
arbitrary  smooth  interfaces,  and  for  a  flat  observation  surface,  the  common  shot 
inversion  algorithm  is  (Dong  [1989],  equation  [3.5.5]) 

B(x)=  \d^K(xsXx)Dm(xs+xr,%),  (1) 

C(x)J 

where 


Vo,  +  a,  Vcos|3(x,  )cosp(^)Vap(^)/a^ 

£(jc-,q,jc)  = - - - - 

T(xs,x  )T  (t  x  )V3pU,)/d? 


(2) 


oo 

=  -j-[Re-Im]f  df  {f e~W  D(f  ,$)  . 
Af  Jo 


(3) 


The  variables  in  these  formulas,  and  those  immediately  following,  are  defined  in  Table 
I.  Except  for  the  factor,  Af,  this  is  just  the  “reflectivity  function”  called  (^(jc)  in 
Bleistein  (1987a).  (We  avoid  the  use  of  p  for  reflectivity  functions  here  because  of  the 
use  of  this  variable  for  incidence  angles  on  the  upper  surface.)  The  introduction  of  the 
scale  Af,  means  that  in  the  neighborhood  of  a  reflector,  the  function  B(x)  has  the 
approximation, 


B(x)  =  R(x,0)yo(x)h(x) . 
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(4) 


1 

Table  I 

R(x,Q) 

angularly  dependent  geometrical  optics  reflection 
coefficient  at  output  point  x . 

YoU) 

band-limited  wavelet,  or  singular  function,  of  unit  height. 

Dlf.Q 

Fourier  transform  of  time  trace  at  receiver  location  q. 

0,(1 ,5) 

Modified  time  section  defined  in  (3). 

cos0 

cosine  of  the  incident  angle  at  reflection  point. 

c(x ) 

background  propagation  speed  at  location  x . 

ray  trace  running  parameter  from  output  point 
to  source  and  receiver,  respectively. 

P(^) 

angles  between  ray  and  upward  vertical  at  xs  and 

ap(^)/a^,  apfx^/a^ 

in-plane  ray  spreading  factors  at  location  ^  and  xs . 

traveltimes  between  output  point  at  x  and  source  at  xs 
or  receiver  at  respectively. 

T (xs,x),  T (^,x) 

product  of  transmission  effects  at  each  interface 
between  output  point  x  and  source  at  xs 
or  receiver  at  %  (  see  reference,  equation  3.5.2). 

V 1 

phase  compensator  due  to  2.5D  (Vif )  is  achieved 
through  this  multiplier  and  the  use  of  Re  -  Im 
combination  of  half  inverse  transform. 

Ai _ 

Area  of  bandpass  filter  applied  to  D  (f ,  £). 

In  this  equation,  h(x )  is  a  slowly  varying  function  of  x,  on  the  scale  of  wavelengths, 
with  /t(x)  =  1  on  the  reflector,  itself.  The  other  new  variables  here  are  also  defined  in 
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Table  I. 


From  these  equations  alone,  one  cannot  estimate  the  change  in  propagation  speed 
across  a  reflector.  It  is  necessary  to  determine  the  distinguished  value  of  0  in  R  (x ,  0) 
to  extract  information  about  the  propagation  speed  of  the  lower  medium.  With  little 
extra  effort,  the  program  also  outputs  the  angularly  dependent  reflection  coefficient 
multiplied  by  cos0  on  the  peak  of  the  singular  function.  The  formula  for  this  second 
inversion  is 


BcU)=~jdkK,(xs,U)Dm(xs+lr.S,)  ,  (5) 

where 

Kc(xs,^,x)  =  K(xs,^,x)c(x)  IV(t,  +tr)l/2  ,  (6) 

with  the  new  variables  again  defined  in  Table  I.  Except  for  the  factor,  c  (or  )/2Aj? ,  this 
is  the  reflectivity  function,  P(jt),  in  Bleistein  (1987a).  In  the  neighborhood  of  a 
reflector. 


Bc(x)  =  R(x,d)  cos  0  Yq(x  )h(x).  (7) 

It  can  now  be  seen  that  the  ratio  of  the  peak  values  of  the  two  outputs, 
Bc  (x  )/B  (x )  =  cos  0,  and  then  either  of  the  outputs  can  be  used  to  estimate  the  change 
in  propagation  speed  across  a  reflector  from  the  dependence  of  the  reflection 
coefficient  on  the  two  propagation  speeds  and  0: 
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cos0 

1 

sin20 

1/2 

c( X) 

c2(x) 

C2U) 

cos0 

1 

sin2© 

ii 

1 

C( X) 

_C+2U) 

c2(x) 

In  this  equation,  c+{x)  is  the  propagation  speed  below  the  reflector  to  be  determined 
from  this  function  and  the  peak  value  of  B  (x ),  for  example,  on  the  reflector. 

In  Bleistein  (1987b),  the  extension  of  this  theory  to  the  variable  density  case  is 
discussed.  If  the  background  density  is  constant,  equations  (1)  and  (5)  still  apply  and 
we  need  only  interpret  the  output  in  terms  of  the  variable  density  geometrical  optics 
reflection  coefficient.  If  the  background  density  on  the  upper  surface  is  not  constant, 
then  one  need  only  modify  those  formulas  by  a  multiplier  of  V p(Xj)/p(x»),  with  p(x) 
the  background  density.  Now  it  is  necessary  to  estimate  the  changes  in  two  medium 
parameters  across  each  reflector.  At  least  two  common  shot  inversions  imaging  the 
same  point  on  a  reflector,  with  different  incidence  angles,  0,  are  necessary  in  this  case 
to  estimate  parameter  changes.  Parsons  (1986)  further  extended  this  scalar  theory  by 
interpreting  the  output  as  a  PP  reflection  coefficient  in  an  elastic  medium.  Now,  at 
least  three  common  shot  inversions  imaging  the  same  point  on  a  reflector  with 
different  incidence  angles,  0,  axe  necessary  to  estimate  parameter  changes.  In  practice, 
all  outputs  imaging  the  same  point  are  used  to  minimize  the  effects  of  noise. 

The  program  CXZCS  implements  equations  (1)  and  (5)  above.  The  filtering 
operation  that  produces  the  modified  data  Dm(x,2;),  described  by  equation  (3),  is 
carried  out  by  a  separate  program.  This  filtering  step  is  independent  of  the  inversion 
process  and  does  not  have  to  he  redone  to  test  different  background  models.  This  step 
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is  carried  out  using  a  fast  Fourier  transform  routine,  producing  Dm(t,%)  on  a  uniform 
grid  in  t.  The  value  Dm(xs  +xr,Z}),  needed  in  equations  (1)  and  (5),  is  then 
approximated  by  three  point  interpolation. 

The  basic  CXZCS  algorithm  (without  interpolation)  is  summarized  below: 


For  each  output  trace 

For  each  depth  z  on  the  output  trat  e 
For  each  receiver  ^ 

T race  rays  from  x  to  q 
Calculate  x's ,  (3'r ,  &s ,  dp's 

Weight  modified  data  at  time  Xs  +  tr  and  sum  into  output  at  x 
Next  receiver 

Multiply  output  at  x  by  constants 
Next  output  depth 
Write  output  trace  to  disk 
Next  output  trace 


Each  output  trace  in  the  depth  section  is  written  as  a  single  record  and,  as  above, 
has  no  header. 


RAYPATH  CALCULATIONS 

The  raypath  calculation  algorithm  in  CXZCS  is  the  same  as  in  Docherty  (1987a). 
The  raypath  between  output  points  and  receivers  on  the  observation  surface  are 
calculated  using  continuation  in  both  interfaces  and  receivers.  That  is,  a  first  ray  path 
trajectory  from  source  to  a  first  output  point  to  receiver  is  determined  for  flat 
interfaces.  Then,  continuation  in  the  interface  shapes  is  used  to  determine  the  ray 
trajectory  for  the  actual  interfaces.  For  each  continuation,  here  and  in  the  following. 
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Newton’s  method  is  used  to  determine  the  new  trajectory.  The  trajectory  from  the 
source  to  the  output  is  determined  this  way,  as  well.  Continuation  in  receiver  position 
is  then  used  to  determine  the  ray  paths  to  all  receivers. 

INTEGRATION  RANGE  SPECIFICATION  AND  INTERPOLATION 

To  image  a  point  on  a  reflector  we  must  insure  that  the  specularly  reflected 
energy  from  that  point  emerges  within  the  range  of  integration  in  equations  (1)  and 
(5).  If  the  dip  of  the  reflector  varies  across  the  section,  it  is  likely  that  the  required 
range  of  integration  will  also  vary.  To  be  safe,  one  usually  specifies  the  range  of  the 
integration  to  include  all  the  receivers  used  in  the  experiment.  This  can  be  costly  and 
may  lead  to  difficulties  and  pathologies  associated  with  ray  tracing. 

For  efficiency,  this  program  allows  the  user  to  define  different  inversion  panels. 
In  each  panel,  the  integration  range  can  be  defined  differently  to  avoid  rays  that  do  not 
give  contributions  to  the  integration.  This  requires  some  a  priori  knowledge  of  the 
true  earth  structure. 

Our  knowledge  of  true  earth  is  expressed  in  the  trial  depth  model  that  we  input  to 
the  program.  Given  the  shot  and  receiver  positions,  this  model  serves  as  our  guide  in 
selecting  the  integration  range.  At  different  output  points,  we  simply  choose  the  limits 
of  integration  (in  terms  of  receiver  numbers)  to  be  used  in  the  integrals  in  equations 
(1)  and  (5).  At  output  points  between  specified  output  locations  the  program  adjusts 
those  limits  so  that  the  integration  range  varies  smoothly  across  the  section.  This  is 
especially  useful  in  regions  of  complex  geological  structure. 
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The  geometry  of  the  seismic  experiment  has  to  be  given  before  the  integration 
range  can  be  specified.  In  CXZCS,  all  receivers  are  located  on  a  flat  observation 
surface,  although  the  theory  allows  for  a  variable  height  upper  surface.  All  the 
receiver  locations  correspond  to  the  trace  locations  of  the  data.  The  user  specifies  the 
number  of  receivers,  the  location  of  the  first  receiver,  and  the  spacing  between  them 
(therefore,  uniform  spread  only).  The  inverted  (migrated)  depth  section  is  a 
rectangular  grid  located  somewhere  within  the  trial  model. 

CXZCS  correctly  treats  geometrical  caustics,  the  bow-tie-like  reflections  that 
appear  on  some  shot  profiles.  In  CXZCS ,  every  possible  ray  is  traced  from  the  output 
point  to  every  receiver.  This  decomposes  the  bow-tie  reflections  into  outputs  on 
different  locations  along  that  interface.  Therefore,  the  inversion  outputs  from  the 
bow-tie-like  reflections  are  still  true  amplitude.  On  the  other  hand,  CXZCS  cannot 
correctly  treat  caustics  of  the  Green’s  function  of  the  background  medium.  In  that 
case,  we  would  have  to  account  for  the  phase  change  in  the  inversion.  This  restriction 
constrains  the  user’s  choice  of  trial  depth  models.  The  theory  can  accommodate  this 
anomaly,  but  this  first  program  was  not  designed  to  deal  with  it. 

There  are  other  considerations  in  the  code  that  are  worth  mentioning  here.  The 
first  is  the  antialiasing  filter.  Spatial  aliasing  of  seismic  data  is  caused  by  insufficient 
sampling.  Receiver  spacing  and  temporal  bandwidth  determine  the  maximum 
emergence  angle  of  rays  at  the  upper  surface  beyond  which  data  will  aliased  to  lower 
transverse  wave  number.  We  define  0a  as  the  critical  angle  of  emergence  at  the  upper 
surface,  beyond  which  information  is  aliased;  0a  is  defined  by 
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(Bleistein,  Cohen  and  Hagin  [1985]).  Here,  /max  is  the  maximum  frequency  in  the 
data,  c  is  the  near  surface  layer  velocity,  and  A£  is  the  receiver  spacing.  For  rays 
whose  emergence  angle  is  greater  than  0a,  we  set  the  amplitude  contribution  of  these 
rays  to  zero,  before  inversion. 

The  second  consideration  is  the  sampling  rate  of  output  depth.  We  take  a  sine 
function  as  a  prototypical  band-limited  delta  function  as  might  appear  in  the  output  of 
the  inversion.  We  then  require  that  the  sampling  rate  in  depth  be  such  that  there  are 
four  sample  points  in  the  interval  of  the  main  lobe  of  a  sine  function  with  the  same 
bandwidth  as  the  given  signal  of  the  filtered  data.  Therefore,  if  A z  is  the  sample 
interval,  one  can  show  that 


Az  = 


_ c _ 

8</m«+/nun) 


(10) 


(See  Dong  [1989]  for  a  discussion.)  We  remark  also  that  the  first  zero  of  the  sine 
function  away  from  its  maximum  belongs  to  the  upper  medium.  That  is,  where  the 
background  medium  changes,  we  process  with  the  upper  medium  propagation  speed 
for  sufficient  depth  to  include  the  main  lobe  of  the  sine  function.  For  deeper  points, 
rays  and  traveltimes  are  calculated  with  the  background  interface  in  place.  The  range 
below  a  discontinuity  of  the  background  for  which  we  use  the  upper  medium 
propagation  speed  is 


c 


(11) 


A*l  = 


2(f  max  +  /  mm) 


The  algorithm  we  have  outlined  so  far  calculates  raypaths  from  every  output  point 
up  to  receivers  within  a  specified  range  on  the  observation  surface.  Unfortunately,  for 
large  data  sets  the  time  spent  on  ray  tracing  can  be  excessive  and  may  account  for  as 
much  as  99%  of  the  total  run  time.  The  interpolation  procedures  we  describe  next 
were  introduced  by  Docherty  (1987a)  in  an  effort  to  reduce  the  cost  of  the  ray  tracing. 
He  found  that  in  most  cases  accurate  images  and  at,  ..utudes  could  be  obtained  from 
only  a  sparse  set  of  ray  path  calculations,  thus  reducing  the  run  time  considerably. 
Typically,  such  calculations  will  be  carried  out  only  for  every  fifth  (or  tenth)  receiver 
point  and  for  each  fifth  (or  tenth)  output  point  both  laterally  and  vertically,  thereby 
saving  a  factor  of  125  to  1000  on  the  ray  path  calculations.  Some  care  is  necessary 
near  the  “top”  of  the  (pseudo)  hyperbola  of  the  traveltime  curve  for  a  particular 
output  point.  Except  for  that,  interpolation  is  equivalent  to  replacing  the  given 
background  model  by  a  nearby  model.  Since  the  background  is  an  approximation,  at 
best,  a  nearby  approximation  would  seem  not  to  be  a  serious  compromise  for  the  CPU 
time  that  is  gained.  For  further  details,  the  reader  is  referred  to  Dong  (1989). 


SYNTHETIC  EXAMPLE 

This  example  demonstrates  the  inversion  on  a  syncline  consisting  of  four 
interfaces.  The  shot  is  located  at  the  middle  of  the  receiver  array.  Sample  ray  paths 
and  the  shot  record  are  displayed  in  Figures  1  and  2,  respectively.  In  the  latter,  we 
can  see  the  bow-tie-like  reflection  from  the  lowest  interface.  The  inversion  result  is 
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shown  in  Figure  3.  We  can  see  that  for  portions  that  rays  have  illuminated,  the 
inversion  result  is  close  to  the  original  model.  Note  that  the  bow-tie  and  attendant 
phase  shifted  source  signature  have  been  successfully  unravelled  by  this  processing. 
All  the  tails  on  the  inversion  section  are  due  to  the  arbitrary  truncation  of  the  data  and 
the  impulse  response  of  the  inversion  operator.  Also,  we  note  that  for  the  portion 
poorly  illuminated  by  rays,  the  inversion  result  is  weak.  Therefore,  we  can  not  expect 
to  determine  the  reflection  coefficient  in  this  part  of  the  output.  Even  in  the  region  of 
illumination,  it  is  not  clear  that  the  output  comes  from  a  sufficient  input  aperture  to 
provide  good  numerical  accurary  (Cohen,  1989).  Estimates  in  regions  of  strongest 
illumination  on  the  syncline  exhibit  errors  up  to  7%  in  the  change  in  propagation 
speed. 


THE  MARATHON  MODEL 

The  data  were  collected  for  Marathon  Oil  Company  at  the  Seismic  Acoustics 
Laboratory  at  the  University  of  Houston  in  a  water  tank  over  a  block  model.  Figure  4 
is  a  cross-section  of  the  model.  The  model  varies  only  slightly  in  the  y-  (out-of¬ 
plane)  direction  so  that  the  2.5D  assumption  used  for  deriving  the  inversion  operators 
in  equations  (1)  and  (5)  are  reasonable  for  these  data.  The  top  layer  was  water.  The 
other  layers  were  various  epoxy  resins.  The  dimensions  in  Figure  4  are  labeled  in 
scaled  feet;  the  original  tank  model  was  only  about  three  feet  wide.  For  brevity,  we 
use  feet  and  seconds  instead  of  scale  feet  and  scale  seconds  in  the  discussion  below. 
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The  data  consist  of  290  shot  records.  There  were  48  receivers  in  an  end-on 


spread.  The  receivers  were  to  the  right  of  the  shot.  The  near  receiver  offset  was  800 
feet  and  the  receiver  spacing  was  80  feet  The  far  receiver  offset  was,  therefore,  4560 
feet.  The  shot  spacing  was  also  80  feet.  The  first  shot  was  located  at  x  =0  feet  and 
the  last  was  at  x  =23200  feet  The  shots  and  receivers  were  at  depth  2=0  feet  shown 
in  Figure  4.  This  is  not  the  water  surface  though.  The  shots  and  receivers  were 
submerged  sufficiently  so  that  no  reflections  from  the  water  surface  were  recorded. 
For  each  shot,  two  seconds  of  data  were  recorded,  sampled  at  4  ms. 

Figure  5  is  a  sample  shot  record  from  these  data.  AGC  has  been  applied  to  this 
record  so  that  more  events  are  visible.  The  shot  location  is  x  =  2000  feet  and  the 
receiver  spread  extends  from  x  =2800  to  x  =6560  feet  The  earliest  event  is  the  direct 
wave.  The  first  curved  event  is  the  water  bottom  reflection.  The  next  two  events  are 
reflections  from  the  second  and  third  interfaces.  The  strong  event  near  1.45  seconds  is 
a  reflection  from  the  model  bottom.  The  reflection  from  the  sawtooth  interface  does 
not  produce  an  easily  identifiable  event  on  this  record. 

The  data  have  been  inverted  three  times  with  different  parameters. 

First  Marathon  inversion 

As  with  migration,  inversion  requires  a  reference  velocity  field.  The  velocities 
used  are  the  exact  velocities  shown  in  Figure  4.  The  interfaces  in  the  derived 
background  model  are  located  nearly  in  the  same  positions  as  the  interfaces  in  Figure 
4.  Differences  occur  where  the  cubic  spline  fit  to  the  control  points  causes  unwanted 
bumps  on  the  interfaces.  The  most  severe  of  these  occurs  on  the  down  thrown  side  of 
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the  fault  cutting  the  third  interface. 


For  each  shot  record,  the  inversion  routine  produced  48  output  depth  traces.  The 
first  output  trace  is  located  at  the  shot  position.  The  trace  separation  is  160  feet.  Each 
trace  has  301  samples.  The  depth  sampling  rate  is  40  feet.  For  each  shot  record,  the 
inversion  produced  an  image  of  the  subsurface  in  the  rectangular  region  from  0  to 
12000  feet  deep  and  from  the  shot  location  to  7520  feet  right  of  the  shot  location. 

The  traveltime  and  amplitude  functions  required  by  the  inversion  formula  are 
computed  by  ray  tracing  from  the  output  point  to  the  source  and  to  the  receivers. 
Rays  were  not,  however,  traced  to  every  output  point  and  every  receiver.  Instead,  rays 
were  traced  from  each  depth  point  on  every  fifth  output  trace  and  to  every  fifth 
receiver.  The  traveltime  and  amplitude  functions  were  interpolated  from  the  values 
obtained  by  the  ray  tracing. 

Figure  6  is  the  inversion  of  the  shot  record  shown  in  Figure  5.  (Note  that  Figure 
5  shows  the  shot  record  after  AGC.  The  inversion  is  applied  to  the  ungained  record.) 
The  inversion  is  a  partial  image  of  the  subsurface.  Shallower  than  1000  feet  is  noise 
from  the  direct  arrival.  The  events  at  3000  feet,  4500  feet,  and  6300  feet  are  images 
of  the  first  three  interfaces.  The  event  at  11200  feet  corresponds  to  the  model  bottom. 
The  fourth  reflector,  the  sawtooth,  is  faint  and  located  at  a  depth  of  9000  feet  and 
distance  of  6000  feet. 

An  inversion  similar  to  Figure  6  is  obtained  for  all  290  shot  records.  Each  yields 
a  different  partial  image  of  the  subsurface.  The  inversions  are  sorted  on  the  output 
trace  location  and  stacked  to  form  a  full  image  of  the  subsurface.  Figure  7  is  a  stack 
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of  all  the  individual  shot  inversions.  Since  the  amplitudes  are  only  significant  in  the 
individual  shot  inversions  and  not  the  stacked  section,  AGC  has  been  applied  so  that 
all  reflections  are  visible.  All  reflectors  are  located  correctly  including  all  teeth  of  the 
fourth  reflector.  There  are  some  short  comings  in  the  inversion. 

The  first  problem  is  that  the  steep  flanks  of  the  dome  are  not  imaged  well.  To 
understand  why,  consider  an  experiment  in  a  constant  velocity  medium  with  a  single 
reflector  and  a  single  source  and  receiver.  The  envelope  of  ail  reflectors  having  the 
same  reflection  time  is  the  familiar  reflection  ellipse  with  the  source  and  receiver  at  the 
foci.  If  the  reflector  has  zero  dip,  the  specular  point  lies  below  the  midpoint  of  the 
source  and  receiver.  As  the  reflector  dip  increases,  the  specular  point  moves  farther  up 
dip  and  laterally  away  from  the  midpoint.  Consequently,  imaging  in  a  region  near  the 
source  and  receiver,  as  with  this  particular  example  inversion,  discriminates  against 
steep  dips. 

Another  problem  with  the  inversion  is  a  phase  reversal  on  the  down  thrown  side 
of  the  fault  cutting  the  third  interface.  There  is  also  a  streak  of  noise  extending  below 
this  phase  reversal.  The  phase  reversal  and  the  noise  beneath  it  are  due  to  the  bump 
on  the  input  model  mentioned  earlier.  The  next  two  inversions  were  performed  to  try 
to  remedy  these  shortcomings. 

Second  Marathon  inversion 

For  the  second  example  inversion  of  the  Marathon  data,  several  parameters  were 
changed.  The  inversion  output  for  each  shot  consisted  of  300  traces  spaced  at  80  feet. 
This  covers  the  entire  model  rather  than  the  limited  portion  of  the  first  inversion. 
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Recorded  energy  from  all  dips  should  be  imageable.  The  background  velocity  model 
was  also  changed;  the  velocities  were  the  same  but  the  interfaces  have  been  changed 
by  reducing  the  number  of  control  points  on  the  interfaces.  Consequently,  the 
interfaces  did  not  have  the  extraneous  bumps  that  the  previous  model  had.  The 
remaining  parameter  changes  concerned  the  traveltime  and  amplitude  function 
interpolation.  Rays  were  traced  from  each  depth  point  on  every  other  output  trace  to 
every  other  receiver.  This  was  done  to  achieve  as  accurate  amplitude  as  possible.  The 
second  inversion  of  all  shot  records  took  3.5  times  as  much  CPU  time  as  the  first 
inversion. 

Inversions  from  all  290  shot  records  were  stacked  to  form  a  complete  image  of 
the  model.  The  image  is  shown  in  Figure  8.  Compared  to  Figure  7,  there  are  some 
improvements  and  also  some  disappointments.  As  expected,  the  steep  flanks  of  the 
dome  are  better  imaged.  The  phase  reversal  in  the  third  interface  was  also  removed. 
There  also  appears  to  be  a  fault  plane  reflection  on  the  third  interface. 

The  main  disappointment  from  this  inversion  is  the  large  degree  of 
migration  smile  noise.  Migration  smiles  occupy  more  area  in  the  shot  inversions  than 
the  reflector  images  do.  This  results  in  many  more  noise  traces  being  stacked  than 
signal  traces.  Although  the  reflector  images  (signal)  are  higher  order  than  the  noise, 
the  large  quantity  of  noise  traces  being  stacked  can  bring  the  noise  to  the  same  level  as 
the  signal.  The  noise  has  nearly  obliterated  the  fifth  tooth  in  the  fourth  (sawtooth) 
interface.  Selective  windowing  of  the  shot  inversions  before  stacking  can  reduce  the 
noise  level  in  the  final  image. 
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Third  Marathon  inversion 


A  third  inversion  with  CXZCS  was  carried  out  on  a  Cray  2  at  the  Minnesota 
Super  Computer  Center.  The  major  objective  of  this  inversion  was  to  see  the  effect  of 
replacing  the  sharp  velocity  contrast  across  the  third  interface  with  two  thin  layers  (one 
wavelength  at  minimum  frequency).  The  result  is  shown  in  Figure  9.  The  salt  flanks 
(in  a  region  in  which  the  background  velocity  was  the  same  as  in  the  previous 
inversions)  and  the  fault  block  in  the  third  interface  are  well  imaged.  Unfortunately, 
there  is  no  improvement  on  the  imaging  of  the  sawtoothed  structure  as  we  had  hoped 
with  this  smoother  background.  Ray  tracing  from  the  neighborhood  of  the  sawtoothed 
structure  below  and  to  the  right  of  the  fault  block  suggests  severe  scattering  of  energy 
from  this  region.  We  believe  that  this  is  what  we  are  seeing  in  the  breakup  of  the 
image  in  this  region. 

The  inversion  operator  delays  introduction  of  a  new  interface  in  the  background 
velocity  until  three  wavelengths  below  that  interface.  We  do  this  because  the  theory 
requires  the  imaging  of  each  reflector  in  the  background  of  the  upper  medium.  To  test 
this,  we  replaced  the  background  with  one  in  which  the  propagation  speeds  were  kept 
constant  below  the  second  interface.  In  this  case,  we  expect  that  the  image  of  the  third 
interface  should  be  as  in  Figure  9,  with  deeper  images  degraded.  The  output  in  Figure 
10  confirms  this  expectation. 

Parameter  determination 

Parameter  estimation  studies  were  carried  out  with  only  limited  success.  There 
were  difficulties  because  the  source  was  directional  and  not  zero  phase.  There  were 
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also  problems  because  common  shot  data  sets  tend  not  to  be  wide  enough  to  give 
accurate  numerical  output;  common  offset  inversion  would  be  better  for  this  purpose. 
This  problem  increases  with  depth.  Furthermore,  the  medium  was  elastic;  an  acoustic 
theory,  even  with  variable  density,  does  not  account  for  energy  lost  to  mode 
conversion,  thus  further  narrowing  the  useful  aperture  of  the  data  for  parameter 
estimation.  The  parameter  that  seemed  to  be  least  affected  by  these  problems  was 
cos0,  since  this  uses  a  ratio  of  inversion  outputs,  in  which  the  errors  shift  the 
numerator  and  die  denominator  in  the  same  direction.  For  details,  see  Emanuel 
(1989a).  Certainly,  more  research  needs  to  be  done  on  this  aspect  of  the  theory. 

CONCLUSIONS 

Common  shot,  c(x,z),  prestack  Kirchhoff  inversion  has  been  demonstrated  to 
image  reflectors.  It  also  provides  a  means  for  parameter  estimation  with  limited 
success  (Parsons,  1986)  for  field  data. 

For  the  Marathon  physical  model  data  set  presented  here,  the  reflectors  were  also 
successfully  imaged.  We  have  shown  that  the  inversion  parameters  may  be  chosen  to 
better  image  the  steep  flanks  of  the  shallow  dome  or  the  deeper  sawteeth.  Both, 
however,  may  be  optimized  if  care  is  taken  to  minimize  stacking  of  noise.  The 
kinematics  of  the  inversion  compare  favorably  to  other  migrations  of  the  data.  More 
research  is  required  on  using  this  method  for  computing  changes  in  medium 
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parameters  from  output  amplitude. 
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FIGURE  CAPTIONS 


FIG.  1:  Syncline  model  and  ray  paths.  Upper  layer  velocity  is  5000  ft/sec;  increment  is 
1000  ft/sec  in  each  layer.  Note:  horizontal  and  vertical  scales  of  plot  are  not 
equal. 


FIG.  2:  The  common  shot  data  generated  for  the  model  of  Figure  1.  Power  gain  was 
used  for  this  display. 


FIG.  3:  The  single  shot  inversion  .esult  of  the  data  in  Figure  2. 


FIG.  4:  Marathon  tank  model  with  given  background  velocities.  These  are  scaled 
variables. 


FIG.  5:  Sample  shot  record  from  the  Marathon  data.  The  shot  location  is  x  =  2000  ft. 
The  receivers  extend  from  x  =  2800  ft.  to  x  =  6560  ft.  AGC  has  been  applied  for 
display. 


FIG.  6:  Shot  inversion  of  shot  record  shown  in  Figure  5.  Reflectors  are  partially  imaged. 


FIG.  7:  Stack  of  inversions  of  all  shot  records.  AGC  has  been  applied. 


FIG.  8:  Stack  of  second  inversions  of  all  shot  records.  Migration  noise  is  greater  than  in 
Figure  12  but  steeper  dips  are  imaged  better.  AGC  has  been  applied  to  the 
stack. 


FIG.  9:  Stack  of  third  inversions  of  all  shot  records. 


FIG.  10:  Stack  of  third  inversions  with  background  velocity  kept  constant  below  the 
second  interface.  Compare  image  of  third  interface  here  with  the  image  in 
Figure  9. 
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